Professor: Marcio Argollo

Contato: Sala A1-02, tel. 2629-5773

Referências básicas:

• Manuais de Fortran 77:
  • http://www.fortran.com/F77_std/rjcnf0001.html
  • http://www.idris.fr/data/cours/lang/fortran/f90/F77.html
  • http://www.stanford.edu/class/me200c/tutorial_77/
• Manuais de Fortran 90
  • http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_fortran90.pdf
  • http://www.cisl.ucar.edu/tcg/consweb/Fortran90/F90Tutorial/tutorial.html
• Fortran 90 para programadores de Fortran 77
  • http://www.nsc.liu.se/~boein/f77to90/
• Manuais de C/C++
  • http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_C.pdf
  • http://imada.sdu.dk/~svalle/courses/dm14-2005/mirror/c/
  • http://www.cs.cf.ac.uk/Dave/C/CE.html
  • http://www.ic.unicamp.br/~cmrubira/aacesta/cpp/cpp15.html
• Livros-texto
  • An Introduction to Computer Simulation Methods: Applications to Physical Systems, H. Gould e J. Tobochnik.
  • Física em computadores, Paulo Murilo Castro de Oliveira e Suzana Maria Moss de Oliveira, São Paulo: Livraria da Física, 2000. (Cópia na biblioteca)

Para uma revisão de comandos básicos em Unix/Linux, utilização do gnuplot e programação em C veja o tutorial apresentado na Semana da Física 2009 "Ferramentas de Computação em Física".

As listas deverão ser entregues por email, compactadas no formato tar.gz e com o nome listaN_seunome.tar.gz

1. Somas não são comutativas no computador! Calcule as somas e com i sendo uma variável ponto flutuante de precisão simples e coloque em um mesmo gráfico as diferencas vs N e vs N para onde .
2. Imagine o empilhamento de dominós de massa unitária e comprimento 2. O primeiro dominó se encontra com o centro de massa na origem, e a cada passo n deslocamos os dominós ja empilhados de para a direita e adicionamos outro dominó na parte inferior da pilha, também com seu centro de massa na origem. Escreva um programa que diz, para um dado , quantos dominós podem ser empilhados até que o centro de massa ultrapasse a largura do primeiro dominó e a pilha tombe. (Veja uma discussão do problema aqui).

Sendo a derivada de uma função definida por e , o método de Euler (com difereças posteriores) corresponde ao truncamento desta expansão em ordem h. Obtemos assim que a derivada no ponto x é aproximada por .

1. Mostre que o método de Euler com diferenças centradas é exato em , onde h é o passo de discretização da derivada. Nesse método escrevemos
2. Calcule a derivada numérica de no ponto x=3 usando o método de Euler simples com passos de discretização e faça um gráfico mostrando a diferença relativa entre a derivada numérica e o valor exato da derivada para variáveis tipo ponto flutuante com precisão simples e dupla. Entre quais valores de h a aproximação é aceitável? (Solução aqui)
3. Mostre que o método de Euler é instável para qualquer passo h na solução do problema de crescimento exponencial
4. Resolva numericamente a equação diferencial com com o método de Euler simples para diferentes valores de passo e compare cada solução nos instantes t=1,2,3,4 e 5 com a solução exata . Voce pode encontrar essa solução exata digitando a equação diferencial no site do WolframAlpha. Experimente!
5. Faça um gráfico comparando a solução numérica com passo e a solução exata entre t=0 e t=5.
6. Resolva a equação acima com o método de Euler com diferenças centradas e os mesmos valores de h. Para obter o primeiro passo use o método de Euler simples com passo de discretização dez vezes menor que o usado no resto do intervalo.

• Referências
  • http://mathworld.wolfram.com/NumericalDifferentiation.html
  • http://davinci.if.ufrgs.br/wiki/index.php/M%C3%A9todo_de_Euler

Ver Reversible multiple time scale molecular dynamics (pdf aqui)

Estes métodos procuram criar esquemas de discretização da evolução Hamiltoniana de um sistema conservativo, de modo que o volume de pontos representativos no espaço de fase seja preservado durante a evolução temporal. Estes esquemas envolvem fatorização do operador Liouvilleano.

1. Use o método de Euler-Cromer e o método de Euler com passo e para resolver numericamente o problema do pêndulo simples na aproximação de pequenas oscilações. Ulitize a condição inicial e . Faça um gráfico contendo as soluções analítica e numérica da a evolução temporal da posição angular por aproximadamente 10 períodos do movimento. Na aproximação numérica use os métodos de Euler e Euler-Cromer. Faça outro gráfico ilustrando a evolução temporal da energia no mesmo alcance temporal para as 3 soluções.

A aproximação para a derivada, como sugerida no método de Euler, introduz erros sistemáticos que acarretam na destruição de certos invariantes da dinâmica original, como a energia total no movimento do pêndulo simples. Para reduzir os efeitos da aproximação numérica à derivada, pode-se recorrer a métodos com aproximações de ordem superior (o método de Euler é de primeira ordem no passo de discretização h).

Ver: http://mathworld.wolfram.com/Runge-KuttaMethod.html
Atenção
Para sistemas de equações diferenciais, ver: http://www.nsc.liu.se/~boein/f77to90/rk.html

1. Resolva numericamente o problema do pêndulo físico utilizando o método de Runge-Kutta de quarta ordem com passo para reproduzir os resultados do artigo "Comportamento crítico no pêndulo simples" do Paulo Murilo. Faça um relatório seguindo os passos deste trabalho. (veja aqui um programa em C e outro em Fortran77 para calcular a evolução temporal do pêndulo)

Nesta sequência de problemas utilize o método de Runge-Kutta de quarta ordem com passo

Considere o pêndulo simples em um meio viscoso, que oferece resistência ao seu movimento com uma força proporcional à sua velocidade . A equação de movimento pode ser escrita como .

1. Faça o gráfico da posição como função do tempo para e com posição angular inicial e . Note que para valores grandes de q devemos recuperar o movimento sem dissipação (conservativo).
2. Considere agora uma força periódica impulsionando o pêndulo com força . A equação de movimento pode ser escrita como .
   1. Construa o gráfico com e e , respectivamente. Observe a mudança de comportamento ao variarmos a amplitude A da força externa. Simule o pêndulo por 100 ciclos da força externa e despreze os primeiros 10 ciclos (quanto vale o período da força externa neste exemplo?).
   2. Faça um gráfico “estroboscópico”, chamado seção de Poincaré, imprimindo os pontos no espaço de fase correspondentes a intervalos de tempo múltiplos do período da força externa .
   3. Obtenha uma série de pontos onde corresponde ao valor da velocidade angular quando pela n-ésima vez durante a evolução temporal do pêndulo. Despreze os valores obtidos para tempos menores que 10 ciclos da força externa. Construa com esses pontos outra seção de Poincaré em um gráfico para os mesmos parâmetros do problema anterior.

Leia aqui sobre o problema.

Considere o mapa logístico

1. Faça os gráficos (escala lin-log) da evolução temporal de x a partir de um valor inicial arbitrário para e .
2. Fazer “gráficos de escada” para os mesmos valores de parâmetro e condição incial x=0.6. Como critério de parada use a condição .
3. Fazer “gráficos de escada” para o caso e condição inicial x=0.3 e para a segunda iterada f(f(x)) com mesma condição inicial.
4. Construir o diagrama de bifurcação do mapa logístico na região .
5. Obter o diagrama de bifurcação para os mapas

a)

b)

ambos na região

1. Obter analiticamente os pontos fixos para órbitas de período 2 e período 4 e analisar as regiões de estabilidade destes pontos fixos. Para encontrar as raízes voces podem usar o WolframAlpha. Escreva, por exemplo, roots of l*x*(1-x)=x
2. O mapa logístico com é ergódico, o que significa que podemos substituir médias temporais por médias sobre configurações de um ensemble e, nesse caso,

onde é a fração de vezes que a órbita passa entre x e x+dx. Calcule, subdividindo o intervalo [0,1] em 100 subintervalos, o histograma para as primeiras 100000 iterações do mapa logístico com e encontre . Esse caso particular admite a solução analítica . Compare os resultados mostrando-os em um mesmo gráfico.

1. Obter os expoentes de Lyapunov para os mapas

a) na região .

b)

c)

ambos b) e c) na região

1. Voce encontra neste arquivo a sequência dos primeiros 10 mil dígitos de e aqui a sequência das primeiras 10 mil bases de nucleotídeos compondo o cromossomo I do genoma do verme Caenorhabditis elegans (C.elegans), o primeiro eucarioto a ter seu genoma completamente sequenciado (ver http://users.rcn.com/jkimball.ma.ultranet/BiologyPages/C/Caen.elegans.html e http://en.wikipedia.org/wiki/Caenorhabditis_elegans para mais detalhes). Neste arquivo 0,1,2,3 correspondem a A,T,C e G, respectivamente.

• Calcular a frequência P(i) de ocorrência de cada símbolo i na sequência.
• Calcular a frequência P(i,j) de ocorrência consecutiva de cada par de símbolos (i,j) na primeira metade da sequência. Utilize essa informação para adivinhar, na segunda metade da sequência, qual o dígito j seguinte a um dado símbolo i da sequência. Qual a sua frequência de acertos? Ler aqui sobre probabilidade condicionada.

1. A função de correlação é uma das ferramentas tradicionais de análise de sinais. Ela mede o grau de correlação entre o valor do sinal medido em um instante e outro sinal medido em um instante de tempo subsequente. Note que, em sinais não estacionários, a função de correlação pode depender do instante inicial da medida . Não trataremos desses casos aqui.
2. Calcule a função de correlação para uma série de 10000 pontos obtidos com a evolução do mapa logístico com (ponto de acumulação, infinitas órbitas periódicas instáveis emergem neste ponto), , (transição caos-ordem em uma janela de período 3) e (ergodicidade comprovada para esse valor de parâmetro). Considere os pontos do mapa logístico com precisão de 5 casas decimais.
3. Calcule a função de correlação para os dígitos de pi e bases de nucleotídeos considerando apenas valores -1 e 1 para digitos abaixo e acima de 5, respectivamente, e para nucleotideos abaixo e acima de 2, respectivamente.

Para ler sobre automata celulares:
http://mathworld.wolfram.com/ElementaryCellularAutomaton.html
http://classes.yale.edu/fractals/CA/welcome.html

1. Desenvolva um código para ilustrar a evolução da regra 90 utilizando somente operações lógicas e um bit para cada sítio do sistema, implementando um sistema de 32 ou 64 bits (dependendo do tamanho da palavra no computador utilizado). Implemente condições de contorno periódicas e livres e condição inicial com um bit central 1 e todos os outros 0.
2. Desenvolva um código para receber o numero de uma dada regra de Wolfram como parâmetro de entrada e retornar a evolução temporal de uma condição inicial com 1 sítio central 1 os outros 0 por T passos de tempo em um sistema com N sítios. Utilize uma palavra de computador por sítio. Mostre graficamente as evoluções temporais para as regras 30, 60, 90, 102, 126 e 150 com N=T=200. Utilize como condição inicial um único sítio (central) com estado s=1 e todos os outros com s=0 e implemente condições de contorno periódicas (s[N]=s[0]) e livres (s[N]=s[-1]=0). Ilustre os resultados graficamente.
3. Dada uma regra de automata (repita para as regras 30, 90, 102 e 126), condições de contorno nula e uma condição inicial com 1 sítio central ocupado, escolha um sítio i e acompanhe sua evolução temporal 0,1,0,0,1,….
   

• Calcule o número médio de bits 1 ao longo da evolução com N=100000 e T=100, 1000, 5000, 10000, 50000 e 100000.
  
• Calcule a função de correlação desta série de bits
  
• Agrupe os bits na evolução em blocos de 8 bits. Transforme o bloco de bits em um inteiro de 8 bits com sinal. Calcule o valor médio do número obtido. Calcule a função de correlação para a série gerada com tal procedimento. Faça um gráfico com a evolução temporal dos números de 8 bits. Para lembrar da representação de inteiros com sinal leia:
  

http://www.cs.cornell.edu/~tomf/notes/cps104/twoscomp.html
http://en.wikipedia.org/wiki/Two's_complement

Números gerados com o computador, uma máquina de Turing determinística, não pode gerar aleatoriedade com seus elementos lógicos. Podemos, no máximo, gerar sequências de números que passem por determinados testes estatísticos de aleatoriedade. A estes números chamamos pseudo-aleatórios e sua obtenção depende crucialmente do aumento de entropia gerado pela perda de bits (informação) nas operações lógicas - para os mais interessados: http://www.springerlink.com/content/jn7x3365386phn46/

Existem geradores de números verdadeiramente aleatórios, baseados na estocasticidade de processos físicos como emissão de fótons em processos radiativos ou espalhamento em divisores de feixes (veja aqui um exemplo). Esses geradores são usados para operações críticas, como criptografia e bingo.

Uma boa fonte de informações a respeito de números pseudo-aleatórios é a biblioteca GSL, que possui, dentre diversas outras aplicações científicas como cálculo de autovalores e solução de sistemas de equações lineares, rotinas que geram números pseudo-aleatórios com diversas receitas. http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Generation.html

1. Usando o gerador congruencial x=(x*A)mod(B) gere 1000 pontos e visualize seu resultado em um espaço 2D considerando pares de numeros aleatórios consecutivos representando pontos no espaço. Utilize A=16807 e 65539 e B=2147483648. Estes números podem ser gerados sem a operação mod com números inteiros de 32 bits com sinal. Visualize seus resultados normalizando o número aleatório para gerar pontos no alcance [-1…1]
2. Usando o mesmo gerador com A=16807 e A=65539, crie pontos em um espaço 3D com cada ponto no espaço obtido a partir de uma trinca de números aleatórios gerados consecutivamente. Visualize seus resultados e, girando a figura no gnuplot, tente visualizar a regularidade nos pontos gerados com 65539. Esses pontos, de acordo com a teoria, caem em 15 planos! (Ver mais aqui.
3. Em 2D, visualize x versus 9x-6y+z onde x,y,z são a trinca de números aleatórios correspondente a um ponto no espaço 3D. Veja como o gerador 65539 é ruim.
4. Gere pares de pontos (x,y) aleatoriamente distribuídos em 2D (com x e y entre -1 e 1). Calcule a fração de números gerados que satisfazem . Essa razão deve ser aproximadamente a razão entre a área de um círculo de raio unitário e um quadrado de lado 2. Esta razão vale . Multiplique sua fração de números por 4 e voce terá calculado pi pelo Método de Monte Carlo. Faça um gráfico mostrando a convergência do valor estimado com o método de Monte Carlo para o valor aproximado de 3.14159265358979323846 (por exemplo, mostre no gráfico a diferença entre o valor aproximado e o valor estimado) utilizando o gerador com A=16807 e compare o resultado obtido com o gerador A=65539. Quantos números aleatórios devemos gerar para se obter de precisão na sua estimativa com A=16807?
5. Ruína do jogador: Considere dois jogadores, o primeiro com n1 moedas e o segundo com n2 moedas em um jogo de cara ou coroa. Cada vez que um jogador perde, ele entrega uma moeda para o outro jogador e o jogo se repete até um dos jogadores perder todo o dinheiro. Realize esse jogo um milhão de vezes e calcule a fração de vezes que cada jogador perdeu todo o dinheiro. Esse valor deve ser interpretado como a probabilidade de cada jogador ir à falência. Compare seu resultado com os valores teóricos para as e . Para ler mais sobre o assunto clique aqui